Let us set some global options for all code chunks in this document.

# Set seed for reproducibility
set.seed(1982) 
# Set global options for all code chunks
knitr::opts_chunk$set(
  # Disable messages printed by R code chunks
  message = TRUE,    
  # Disable warnings printed by R code chunks
  warning = TRUE,    
  # Show R code within code chunks in output
  echo = TRUE,        
  # Include both R code and its results in output
  include = TRUE,     
  # Evaluate R code chunks
  eval = TRUE,       
  # Enable caching of R code chunks for faster rendering
  cache = FALSE,      
  # Align figures in the center of the output
  fig.align = "center",
  # Enable retina display for high-resolution figures
  retina = 2,
  # Show errors in the output instead of stopping rendering
  error = TRUE,
  # Do not collapse code and output into a single block
  collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
  fig_count <<- fig_count + 1
  paste0("Figure ", fig_count, ": ", caption)
}
# Define the function to truncate a number to two decimal places
truncate_to_two <- function(x) {
  floor(x * 100) / 100
}

m1table <- rSPDE:::m1table
m2table <- rSPDE:::m2table
m3table <- rSPDE:::m3table
m4table <- rSPDE:::m4table
# install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/testing"), dep=TRUE)
# inla.upgrade(testing = TRUE)
# remotes::install_github("inlabru-org/inlabru", ref = "devel")
# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(INLA)
## Loading required package: Matrix
## This is INLA_25.05.01 built 2025-05-01 18:43:33 UTC.
##  - See www.r-inla.org/contact-us for how to get help.
##  - List available models/likelihoods/etc with inla.list.models()
##  - Use inla.doc(<NAME>) to access documentation
##  - Consider upgrading R-INLA to testing[25.05.07] or stable[24.12.11].
library(inlabru)
## Loading required package: fmesher
library(rSPDE)
## This is rSPDE 2.5.1
## - See https://davidbolin.github.io/rSPDE for vignettes and manuals.
library(MetricGraph)
## This is MetricGraph 1.4.1
## - See https://davidbolin.github.io/MetricGraph for vignettes and manuals.
## 
## Attaching package: 'MetricGraph'
## The following object is masked from 'package:stats':
## 
##     filter
library(grateful)

library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

Utilitary functions

# For each m and beta, this function returns c_m/b_{m+1} and the roots of rb and rc
my.get.roots <- function(order, beta) {
  mt <- get(paste0("m", order, "table"))
  rb <- rep(0, order + 1)
  rc <- rep(0, order)
  if(order == 1) {
      rc = approx(mt$beta, mt[[paste0("rc")]], beta)$y
    } else {
      rc = sapply(1:order, function(i) {
        approx(mt$beta, mt[[paste0("rc.", i)]], beta)$y
      })
    }
    rb = sapply(1:(order+1), function(i) {
      approx(mt$beta, mt[[paste0("rb.", i)]], xout = beta)$y
    })
    factor = approx(mt$beta, mt$factor, xout = beta)$y
  return(list(rb = rb, rc = rc, factor = factor))
}

# Function the polynomial coefficients in increasing order like a+bx+cx^2+...
poly.from.roots <- function(roots) {
  
  coef <- 1
  for (r in roots) {coef <- convolve(coef, c(1, -r), type = "open")}
  return(coef)
}
# Function to compute the partial fraction parameters
compute.partial.fraction.param <- function(factor, pr_roots, pl_roots, cte) {
  
  pr_coef <- c(0, poly.from.roots(pr_roots)) 
  pl_coef <- poly.from.roots(pl_roots) 
  factor_pr_coef <- pr_coef
  pr_plus_pl_coef <- factor_pr_coef + cte/factor * pl_coef
  res <- gsignal::residue(factor_pr_coef, pr_plus_pl_coef)
  return(list(r = res$r, p = res$p, k = res$k)) 
}
# Function to compute the fractional operator
my.fractional.operators.frac <- function(L, beta, C, scale.factor, m = 1, time_step) {
  
  C <- Matrix::Diagonal(dim(C)[1], rowSums(C)) 
  Ci <- Matrix::Diagonal(dim(C)[1], 1 / rowSums(C)) 
  I <- Matrix::Diagonal(dim(C)[1])
  L <- L / scale.factor 
  LCi <- L %*% Ci
  if(beta == 1){
    return(list(Ci = Ci, C = C, LCi = LCi, L = L, m = m, beta = beta, LHS = C + time_step * L))
  } else {
    roots <- my.get.roots(m, beta)
    poles_rs_k <- compute.partial.fraction.param(roots$factor, roots$rc, roots$rb, time_step)

    partial_fraction_terms <- list()
    for (i in 1:(m+1)) {partial_fraction_terms[[i]] <- (LCi - poles_rs_k$p[i] * I)/poles_rs_k$r[i]}
    partial_fraction_terms[[m+2]] <- ifelse(is.null(poles_rs_k$k), 0, poles_rs_k$k) * I
    return(list(Ci = Ci, C = C, LCi = LCi, L = L, m = m, beta = beta, partial_fraction_terms = partial_fraction_terms))
  }
}
# Function to solve the iteration
my.solver.frac <- function(obj, v){
  beta <- obj$beta
  m <- obj$m
  C <- obj$C
  Ci <- obj$Ci
  if (beta == 1){
    return(solve(obj$LHS, v))
  } else {
    partial_fraction_terms <- obj$partial_fraction_terms
    output <- partial_fraction_terms[[m+2]] %*% v
    for (i in 1:(m+1)) {output <- output + solve(partial_fraction_terms[[i]], v)}
    return(Ci %*% output)
  }
}
# Function to build a tadpole graph and create a mesh
gets.graph.tadpole <- function(h){
  edge1 <- rbind(c(0,0),c(1,0))
  theta <- seq(from=-pi,to=pi,length.out = 10000)
  edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
  edges <- list(edge1, edge2)
  graph <- metric_graph$new(edges = edges)
  graph$set_manual_edge_lengths(edge_lengths = c(1,2))
  graph$build_mesh(h = h)
  return(graph)
}
# Function to compute the eigenfunctions 
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2]) 
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2]) 

if(k==0){ 
  f.e1 <- rep(1,length(x1)) 
  f.e2 <- rep(1,length(x2)) 
  f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
  f = list(phi=f1/sqrt(3)) 
  
} else {
  f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2) 
  f.e2 <- sin(pi*k*x2/2)                  
  
  f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1]) 
  
  if((k %% 2)==1){ 
    f = list(phi=f1/sqrt(3)) 
  } else { 
    f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
    f.e2 <- cos(pi*k*x2/2)
    f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1]) 
    f <- list(phi=f1,psi=f2/sqrt(3/2))
  }
}
return(f)
}
# Function to order the vertices for plotting
plotting.order <- function(v, graph){
  edge_number <- graph$mesh$VtE[, 1]
  pos <- sum(edge_number == 1)+1
  return(c(v[1], v[3:pos], v[2], v[(pos+1):length(v)], v[2]))
}
# Function to plot in 3D
graph.plotter.3d <- function(graph, 
                             U_true, 
                             U_approx1, 
                             #U_approx2, 
                             time_seq, 
                             time_step){
  x <- graph$mesh$V[, 1]
  y <- graph$mesh$V[, 2]
  x <- plotting.order(x, graph)
  y <- plotting.order(y, graph)
  weights <- graph$mesh$weights
  
  cumsum1 <- sqrt(time_step * cumsum(t(weights) %*% (U_true - U_approx1)^2))
  #cumsum2 <- sqrt(time_step * cumsum(t(weights) %*%(U_true - U_approx2)^2))
  #cumsum3 <- sqrt(time_step * cumsum(t(weights) %*% (U_approx1 - U_approx2)^2))
  
  U_true <- apply(U_true, 2, plotting.order, graph = graph)
  U_approx1 <- apply(U_approx1, 2, plotting.order, graph = graph)
  #U_approx2 <- apply(U_approx2, 2, plotting.order, graph = graph)
  
  plot_data <- data.frame(
  x = rep(x, times = ncol(U_true)),
  y = rep(y, times = ncol(U_true)),
  z_true = as.vector(U_true),
  z_approx1 = as.vector(U_approx1),
  #z_approx2 = as.vector(U_approx2),
  frame = rep(time_seq, each = length(x)))
  
  x_range <- range(x)
  y_range <- range(y)
  z_range <- range(c(U_true, 
                     #U_approx2, 
                     U_approx1))
  
  p <- plot_ly(plot_data, frame = ~frame) %>%
  add_trace(x = ~x, y = ~y, z = ~z_true, type = "scatter3d", mode = "lines", name = "T",
            line = list(color = "green", width = 2)) %>%
  add_trace(x = ~x, y = ~y, z = ~z_approx1, type = "scatter3d", mode = "lines", name = "A1",
            line = list(color = "red", width = 2)) %>%
  # add_trace(x = ~x, y = ~y, z = ~z_approx2, type = "scatter3d", mode = "lines", name = "A2", 
  #           line = list(color = "blue", width = 2)) %>%
  layout(
    scene = list(
      xaxis = list(title = "x", range = x_range),
      yaxis = list(title = "y", range = y_range),
      zaxis = list(title = "z", range = z_range),
      aspectratio = list(x = 2.4, y = 1.2, z = 1.2),
      camera = list(eye = list(x = 1.5, y = 1.5, z = 1), center = list(x = 0, y = 0, z = 0))),
    updatemenus = list(list(type = "buttons", showactive = FALSE,
                            buttons = list(
          list(label = "Play", method = "animate",
               args = list(NULL, list(frame = list(duration = 100, redraw = TRUE), fromcurrent = TRUE))),
          list(label = "Pause", method = "animate",
               args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
        )
      )
    ),
    title = "Time: 0"
  ) %>% plotly_build()
  for (i in seq_along(p$x$frames)) {
    t <- time_seq[i]
    err1 <- signif(cumsum1[i], 6)
    #err2 <- signif(cumsum2[i], 4)
    #err3 <- signif(cumsum3[i], 4)
    p$x$frames[[i]]$layout <- list(title = paste0("Time: ", t, 
                                                  # " | cum E=T-A2: ", err2, 
                                                  # " | cum E=A1-A2: ", err3,
                                                  " | cum E=T-A1: ", err1))
  }
  return(p)
}
# Function to plot the error at each time step
error.at.each.time.plotter <- function(graph, U_true, U_approx1, U_approx2, time_seq, time_step) {
  weights <- graph$mesh$weights
  error_at_each_time1 <- t(weights) %*% (U_true - U_approx1)^2
  error_at_each_time2 <- t(weights) %*% (U_true - U_approx2)^2
  error_between_both_approx <- t(weights) %*% (U_approx1 - U_approx2)^2
  
  error1 <- sqrt(as.double(t(weights) %*% (U_true - U_approx1)^2 %*% rep(time_step, ncol(U_true))))
  error2 <- sqrt(as.double(t(weights) %*% (U_true - U_approx2)^2 %*% rep(time_step, ncol(U_true))))
  errorb <- sqrt(as.double(t(weights) %*% (U_approx1 - U_approx2)^2 %*% rep(time_step, ncol(U_true))))
  
  fig <- plot_ly() %>% 
  add_trace(
  x = ~time_seq, y = ~error_at_each_time1, type = 'scatter', mode = 'lines+markers',
  line = list(color = 'red', width = 2),
  marker = list(color = 'red', size = 4),
  name = paste0("E=T-A1: ", sprintf("%.3e", error1))
) %>% 
  add_trace(
  x = ~time_seq, y = ~error_at_each_time2, type = 'scatter', mode = 'lines+markers',
  line = list(color = 'blue', width = 2, dash = "dot"),
  marker = list(color = 'blue', size = 4),
  name = paste0("E=T-A2: ", sprintf("%.3e", error2))
) %>% 
  add_trace(
  x = ~time_seq, y = ~error_between_both_approx, type = 'scatter', mode = 'lines+markers',
  line = list(color = 'orange', width = 2),
  marker = list(color = 'orange', size = 4),
  name = paste0("E=A1-A2: ", sprintf("%.3e", errorb))
) %>% 
  layout(
  title = "Error at Each Time Step",
  xaxis = list(title = "Error at Each Time Step"),
  yaxis = list(title = "Error"),
  legend = list(x = 0.1, y = 0.9)
)
  return(fig)
}
# Function to compute the eigenvalues and eigenfunctions
gets.eigen.params <- function(N_finite = 4, kappa = 1, alpha = 0.5, graph){
  EIGENVAL_ALPHA <- NULL
  EIGENFUN <- NULL
  for (j in 0:N_finite) {
      lambda_j_alpha <- (kappa^2 + (j*pi/2)^2)^(alpha/2)
      e_j <- tadpole.eig(j,graph)$phi
      EIGENVAL_ALPHA <- c(EIGENVAL_ALPHA, lambda_j_alpha)         
      EIGENFUN <- cbind(EIGENFUN, e_j)
      if (j>0 && (j %% 2 == 0)) {
        lambda_j_alpha <- (kappa^2 + (j*pi/2)^2)^(alpha/2)
        e_j <- tadpole.eig(j,graph)$psi
        EIGENVAL_ALPHA <- c(EIGENVAL_ALPHA, lambda_j_alpha)         
        EIGENFUN <- cbind(EIGENFUN, e_j)
      }
  }
  return(list(EIGENVAL_ALPHA = EIGENVAL_ALPHA,
              EIGENFUN = EIGENFUN))
}

We want to solve the fractional diffusion equation \[\begin{equation} \label{eq:maineq} \partial_t u+(\kappa^2-\Delta_\Gamma)^{\alpha/2} u=f \text { on } \Gamma \times(0, T), \quad u(0)=u_0 \text { on } \Gamma, \end{equation}\] where \(u\) satisfies the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \left\{\phi\in C(\Gamma)\;\Big|\; \forall v\in V: \sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\} \end{equation}\] The solution is given by \[\begin{equation} \label{eq:sol_reprentation} u(s,t) = \displaystyle\sum_{j\in\mathbb{N}}e^{-\lambda^{\alpha/2}_jt}\left(u_0, e_j\right)_{L_2(\Gamma)}e_j(s) + \int_0^t \displaystyle\sum_{j\in\mathbb{N}}e^{-\lambda^{\alpha/2}_j(t-r)}\left(f(\cdot, r), e_j\right)_{L_2(\Gamma)}e_j(s)dr. \end{equation}\]

If we choose \(w_j\) and \(v_j\) and take the initial condition and the right hand side funciton as

\[\begin{equation} u_0(s) = \sum_{j=0}^{N} w_j e_j(s) \text{ and so } \left(u_0, e_j\right)_{L_2(\Gamma)} = w_j, \end{equation}\] \[\begin{equation} \text{In matrix notation: } \quad\boldsymbol{U}_0 = \boldsymbol{E}_h\boldsymbol{c}, \quad \boldsymbol{E}_h = \left[e_0, e_1, \ldots, e_{N}\right], \quad \boldsymbol{w} = \left[w_0, w_1, \ldots, w_{N}\right]^\top, \end{equation}\] \[\begin{equation} \text{In } \texttt{R}: \texttt{U_0 <- EIGENFUN %*% coeff_U_0} \end{equation}\] and \[\begin{equation} f(s,t) = \sum_{j=0}^{N} v_j e^{-\lambda^{\alpha/2}_jt} e_j(s) \text{ and so } \left(f(\cdot,r), e_j\right)_{L_2(\Gamma)} = v_j e^{-\lambda^{\alpha/2}_jr}, \end{equation}\] \[\begin{equation} \text{In matrix notation: } \quad\boldsymbol{f} = \boldsymbol{E}_h \boldsymbol{V}, \quad \boldsymbol{V}_{ji} = v_je^{-\lambda^{\alpha/2}_jt_i} \end{equation}\] \[\begin{equation} \text{In } \texttt{R}: \texttt{FF <- EIGENFUN %*% (coeff_FF * exp(-outer(EIGENVAL_ALPHA, time_seq)))} \end{equation}\] then the solution is given by \[\begin{align} u(s,t) &= \displaystyle\sum_{j=0}^{N}w_je^{-\lambda^{\alpha/2}_jt}e_j(s) + \int_0^t \displaystyle\sum_{j=0}^Ne^{-\lambda^{\alpha/2}_j(t-r)}v_j e^{-\lambda^{\alpha/2}_jr}e_j(s)dr\\ &= \displaystyle\sum_{j=0}^{N}w_je^{-\lambda^{\alpha/2}_jt}e_j(s) + t \displaystyle\sum_{j=0}^Nv_j e^{-\lambda^{\alpha/2}_jt}e_j(s)\\ &= \displaystyle\sum_{j=0}^{N}w_je^{-\lambda^{\alpha/2}_jt}e_j(s) + tf(s,t)\\ &= \displaystyle\sum_{j=0}^{N}(w_j+tv_j)e^{-\lambda^{\alpha/2}_jt}e_j(s) \end{align}\]

\[\begin{equation} \text{In matrix notation: } \quad\boldsymbol{U} =\boldsymbol{E}_h \boldsymbol{W} + (\boldsymbol{t}\boldsymbol{f}^\top)^\top, \quad \boldsymbol{W}_{ji} = w_je^{-\lambda^{\alpha/2}_jt_i},\quad \boldsymbol{t} = \left[t_0, t_1, \ldots, t_{K}\right]^\top \end{equation}\] \[\begin{equation} \text{In } \texttt{R}: \texttt{U_true <- EIGENFUN %*% (coeff_U_0 * exp(-outer(EIGENVAL_ALPHA, time_seq)))+ t(time_seq * t(FF))} \end{equation}\]

kappa <- 1
alpha <- 0.5 # from 0.5 to 2
m = 1
beta <- alpha/2
N_finite = 4 # choose even
adjusted_N_finite <- N_finite + N_finite/2 + 1
# Coefficients for u_0 and f
coeff_U_0 <- 50*(1:adjusted_N_finite)^-1
coeff_U_0[-5] <- 0
coeff_FF <- rep(0, adjusted_N_finite)
coeff_FF[7] <- 10
T_final <- 2

overkill_time_power <- 9
overkill_h_power <- 9

time_steps <- 0.1 * 2^-c(overkill_time_power:0)
hs <- 0.1 * 2^-c(overkill_h_power:0)

overkill_time_step <- time_steps[1]
overkill_h <- hs[1]

overkill_time_seq <- seq(0, T_final, by = overkill_time_step)

overkill_graph <- gets.graph.tadpole(h = overkill_h)
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.16 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
overkill_graph$compute_fem()
overkill_weights <- overkill_graph$mesh$weights
overkill_eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = overkill_graph)
EIGENVAL_ALPHA <- overkill_eigen_params$EIGENVAL_ALPHA
overkill_EIGENFUN <- overkill_eigen_params$EIGENFUN


overkill_U_true <- overkill_EIGENFUN %*% 
  outer(1:length(coeff_U_0), 1:length(overkill_time_seq), 
        function(i, j) (coeff_U_0[i] + coeff_FF[i] * overkill_time_seq[j]) * exp(-EIGENVAL_ALPHA[i] * overkill_time_seq[j]))

by_vector <- 2^c(0:overkill_time_power)

errors <- matrix(NA, nrow = length(hs), ncol = length(time_steps))
for (i in 1:length(hs)) {
  h <- hs[i]
  graph <- gets.graph.tadpole(h = h)
  graph$compute_fem()
  G <- graph$mesh$G
  C <- graph$mesh$C
  L <- kappa^2*C + G
  eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = graph)
  EIGENFUN <- eigen_params$EIGENFUN
  U_0 <- EIGENFUN %*% coeff_U_0
  A <- graph$fem_basis(overkill_graph$get_mesh_locations())
  for (j in 1:length(time_steps)) {
    time_step <- time_steps[j]
    print(by_vector[j])
    coarse_indices <- seq(1, length(overkill_time_seq), by = by_vector[j])
    time_seq <- overkill_time_seq[coarse_indices]
    my_op_frac <- my.fractional.operators.frac(L, beta, C, scale.factor = kappa^2, m = m, time_step)
    INT_BASIS_EIGEN <- t(overkill_EIGENFUN) %*% overkill_graph$mesh$C %*% A
    FF_approx <- t(INT_BASIS_EIGEN) %*% 
      outer(1:length(coeff_FF), 1:length(time_seq), 
        function(i, j) coeff_FF[i] * exp(-EIGENVAL_ALPHA[i] * time_seq[j]))
    U_approx <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
    U_approx[, 1] <- U_0
    for (k in 1:(length(time_seq) - 1)) {
      U_approx[, k + 1] <- as.matrix(my.solver.frac(my_op_frac, my_op_frac$C %*% U_approx[, k] + time_step * FF_approx[, k + 1]))
      }
    sliced_overkill_U_true <- overkill_U_true[, coarse_indices]
    projected_U_approx <- A %*% U_approx
    errors[i,j] <- sqrt(as.double(t(overkill_weights) %*% (sliced_overkill_U_true - projected_U_approx)^2 %*% rep(time_step, length(time_seq))))
  }
}
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.15 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.37 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.30 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.30 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.51 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.18 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.25 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.35 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.29 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
## Starting graph creation...
## LongLat is set to FALSE
## Creating edges...
## Setting edge weights...
## Computing bounding box...
## Setting up edges
## Merging close vertices
## Total construction time: 0.16 secs
## Creating and updating vertices...
## Storing the initial graph...
## Computing the relative positions of the edges...
## [1] 1
## [1] 2
## [1] 4
## [1] 8
## [1] 16
## [1] 32
## [1] 64
## [1] 128
## [1] 256
## [1] 512
print(errors)
##             [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
##  [1,] 0.04698377 0.04654427 0.04569847 0.04414803 0.04167862 0.03971968
##  [2,] 0.04698315 0.04654364 0.04569783 0.04414737 0.04167791 0.03971891
##  [3,] 0.04698068 0.04654115 0.04569528 0.04414471 0.04167506 0.03971583
##  [4,] 0.04697082 0.04653117 0.04568509 0.04413409 0.04166366 0.03970352
##  [5,] 0.04693144 0.04649136 0.04564441 0.04409171 0.04161815 0.03965433
##  [6,] 0.04677529 0.04633349 0.04548309 0.04392359 0.04143755 0.03945901
##  [7,] 0.04617306 0.04572443 0.04486042 0.04327407 0.04073844 0.03870085
##  [8,] 0.04414494 0.04367132 0.04275719 0.04107114 0.03834809 0.03607640
##  [9,] 0.04316569 0.04266440 0.04169331 0.03988714 0.03690499 0.03412564
## [10,] 0.10767692 0.10745260 0.10702105 0.10622713 0.10492180 0.10347597
##             [,7]       [,8]      [,9]     [,10]
##  [1,] 0.04794803 0.08813812 0.1838800 0.3732394
##  [2,] 0.04794735 0.08813768 0.1838797 0.3732391
##  [3,] 0.04794462 0.08813595 0.1838785 0.3732380
##  [4,] 0.04793371 0.08812900 0.1838738 0.3732337
##  [5,] 0.04789015 0.08810129 0.1838550 0.3732164
##  [6,] 0.04771720 0.08799126 0.1837803 0.3731477
##  [7,] 0.04704621 0.08756444 0.1834883 0.3728760
##  [8,] 0.04471973 0.08607459 0.1824294 0.3718457
##  [9,] 0.04236829 0.08382289 0.1799775 0.3686312
## [10,] 0.10530228 0.12515923 0.1978960 0.3703496
graph.plotter.3d(overkill_graph, sliced_overkill_U_true, projected_U_approx, time_seq, time_step)
library(ggplot2)
library(reshape2)
library(plotly)


# Create matrix
mat <- errors

# Convert to data frame
df <- melt(mat)
colnames(df) <- c("Y", "X", "Value")

# ggplot
p <- ggplot(df, aes(x = X, y = Y, fill = Value, text = paste("Value:", round(Value, 2)))) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "red") +
  coord_fixed() +
  theme_minimal()

p

# Convert to interactive plot
ggplotly(p, tooltip = "text")
image(1:10, 1:10, t(mat)[, nrow(mat):1], col = heat.colors(100), xlab = "", ylab = "")

LS0tCnRpdGxlOiAiU29sdmluZyBhIHBhcmFib2xpYyBlcXVhdGlvbiIKZGF0ZTogIkNyZWF0ZWQ6IDIwLTA0LTIwMjUuIExhc3QgbW9kaWZpZWQ6IGByIGZvcm1hdChTeXMudGltZSgpLCAnJWQtJW0tJVkuJylgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIG1hdGhqYXg6ICJodHRwczovL2Nkbi5qc2RlbGl2ci5uZXQvbnBtL21hdGhqYXhAMy9lczUvdGV4LW1tbC1jaHRtbC5qcyIKICAgIGhpZ2hsaWdodDogcHlnbWVudHMKICAgIHRoZW1lOiBmbGF0bHkKICAgIGNvZGVfZm9sZGluZzogaGlkZSAjIGNsYXNzLnNvdXJjZSA9ICJmb2xkLWhpZGUiIHRvIGhpZGUgY29kZSBhbmQgYWRkIGEgYnV0dG9uIHRvIHNob3cgaXQKICAgIGRmX3ByaW50OiBwYWdlZAogICAgIyB0b2M6IHRydWUKICAgICMgdG9jX2Zsb2F0OgogICAgIyAgIGNvbGxhcHNlZDogdHJ1ZQogICAgIyAgIHNtb290aF9zY3JvbGw6IHRydWUKICAgIG51bWJlcl9zZWN0aW9uczogZmFsc2UKICAgIGZpZ19jYXB0aW9uOiB0cnVlCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlCmFsd2F5c19hbGxvd19odG1sOiB0cnVlCmJpYmxpb2dyYXBoeTogCiAgLSByZWZlcmVuY2VzLmJpYgogIC0gZ3JhdGVmdWwtcmVmcy5iaWIKaGVhZGVyLWluY2x1ZGVzOgogIC0gXG5ld2NvbW1hbmR7XGFyfXtcbWF0aGJie1J9fQogIC0gXG5ld2NvbW1hbmR7XGxsYXZ9WzFde1xsZWZ0XHsjMVxyaWdodFx9fQogIC0gXG5ld2NvbW1hbmR7XHBhcmV9WzFde1xsZWZ0KCMxXHJpZ2h0KX0KICAtIFxuZXdjb21tYW5ke1xOY2FsfXtcbWF0aGNhbHtOfX0KICAtIFxuZXdjb21tYW5ke1xWY2FsfXtcbWF0aGNhbHtWfX0KICAtIFxuZXdjb21tYW5ke1xFY2FsfXtcbWF0aGNhbHtFfX0KICAtIFxuZXdjb21tYW5ke1xXY2FsfXtcbWF0aGNhbHtXfX0KLS0tCgpgYGB7ciB4YXJpbmdhbkV4dHJhLWNsaXBib2FyZCwgZWNobyA9IEZBTFNFfQpodG1sdG9vbHM6OnRhZ0xpc3QoCiAgeGFyaW5nYW5FeHRyYTo6dXNlX2NsaXBib2FyZCgKICAgIGJ1dHRvbl90ZXh0ID0gIjxpIGNsYXNzPVwiZmEtc29saWQgZmEtY2xpcGJvYXJkXCIgc3R5bGU9XCJjb2xvcjogIzAwMDA4QlwiPjwvaT4iLAogICAgc3VjY2Vzc190ZXh0ID0gIjxpIGNsYXNzPVwiZmEgZmEtY2hlY2tcIiBzdHlsZT1cImNvbG9yOiAjOTBCRTZEXCI+PC9pPiIsCiAgICBlcnJvcl90ZXh0ID0gIjxpIGNsYXNzPVwiZmEgZmEtdGltZXMtY2lyY2xlXCIgc3R5bGU9XCJjb2xvcjogI0Y5NDE0NFwiPjwvaT4iCiAgKSwKICBybWFya2Rvd246Omh0bWxfZGVwZW5kZW5jeV9mb250X2F3ZXNvbWUoKQopCmBgYAoKCmBgYHtjc3MsIGVjaG8gPSBGQUxTRX0KYm9keSAubWFpbi1jb250YWluZXIgewogIG1heC13aWR0aDogMTAwJSAhaW1wb3J0YW50OwogIHdpZHRoOiAxMDAlICFpbXBvcnRhbnQ7Cn0KYm9keSB7CiAgbWF4LXdpZHRoOiAxMDAlICFpbXBvcnRhbnQ7Cn0KCmJvZHksIHRkIHsKICAgZm9udC1zaXplOiAxNnB4Owp9CmNvZGUucnsKICBmb250LXNpemU6IDE0cHg7Cn0KcHJlIHsKICBmb250LXNpemU6IDE0cHgKfQouY3VzdG9tLWJveCB7CiAgYmFja2dyb3VuZC1jb2xvcjogI2Y1ZjdmYTsgLyogTGlnaHQgZ3JleS1ibHVlIGJhY2tncm91bmQgKi8KICBib3JkZXItY29sb3I6ICNlMWU4ZWQ7IC8qIExpZ2h0IGJvcmRlciBjb2xvciAqLwogIGNvbG9yOiAjMmMzZTUwOyAvKiBEYXJrIHRleHQgY29sb3IgKi8KICBwYWRkaW5nOiAxNXB4OyAvKiBQYWRkaW5nIGluc2lkZSB0aGUgYm94ICovCiAgYm9yZGVyLXJhZGl1czogNXB4OyAvKiBSb3VuZGVkIGNvcm5lcnMgKi8KICBtYXJnaW4tYm90dG9tOiAyMHB4OyAvKiBTcGFjaW5nIGJlbG93IHRoZSBib3ggKi8KfQouY2FwdGlvbiB7CiAgbWFyZ2luOiBhdXRvOwogIHRleHQtYWxpZ246IGNlbnRlcjsKICBtYXJnaW4tYm90dG9tOiAyMHB4OyAvKiBTcGFjaW5nIGJlbG93IHRoZSBib3ggKi8KfQpgYGAKCgpMZXQgdXMgc2V0IHNvbWUgZ2xvYmFsIG9wdGlvbnMgZm9yIGFsbCBjb2RlIGNodW5rcyBpbiB0aGlzIGRvY3VtZW50LgoKCmBgYHtyfQojIFNldCBzZWVkIGZvciByZXByb2R1Y2liaWxpdHkKc2V0LnNlZWQoMTk4MikgCiMgU2V0IGdsb2JhbCBvcHRpb25zIGZvciBhbGwgY29kZSBjaHVua3MKa25pdHI6Om9wdHNfY2h1bmskc2V0KAogICMgRGlzYWJsZSBtZXNzYWdlcyBwcmludGVkIGJ5IFIgY29kZSBjaHVua3MKICBtZXNzYWdlID0gVFJVRSwgICAgCiAgIyBEaXNhYmxlIHdhcm5pbmdzIHByaW50ZWQgYnkgUiBjb2RlIGNodW5rcwogIHdhcm5pbmcgPSBUUlVFLCAgICAKICAjIFNob3cgUiBjb2RlIHdpdGhpbiBjb2RlIGNodW5rcyBpbiBvdXRwdXQKICBlY2hvID0gVFJVRSwgICAgICAgIAogICMgSW5jbHVkZSBib3RoIFIgY29kZSBhbmQgaXRzIHJlc3VsdHMgaW4gb3V0cHV0CiAgaW5jbHVkZSA9IFRSVUUsICAgICAKICAjIEV2YWx1YXRlIFIgY29kZSBjaHVua3MKICBldmFsID0gVFJVRSwgICAgICAgCiAgIyBFbmFibGUgY2FjaGluZyBvZiBSIGNvZGUgY2h1bmtzIGZvciBmYXN0ZXIgcmVuZGVyaW5nCiAgY2FjaGUgPSBGQUxTRSwgICAgICAKICAjIEFsaWduIGZpZ3VyZXMgaW4gdGhlIGNlbnRlciBvZiB0aGUgb3V0cHV0CiAgZmlnLmFsaWduID0gImNlbnRlciIsCiAgIyBFbmFibGUgcmV0aW5hIGRpc3BsYXkgZm9yIGhpZ2gtcmVzb2x1dGlvbiBmaWd1cmVzCiAgcmV0aW5hID0gMiwKICAjIFNob3cgZXJyb3JzIGluIHRoZSBvdXRwdXQgaW5zdGVhZCBvZiBzdG9wcGluZyByZW5kZXJpbmcKICBlcnJvciA9IFRSVUUsCiAgIyBEbyBub3QgY29sbGFwc2UgY29kZSBhbmQgb3V0cHV0IGludG8gYSBzaW5nbGUgYmxvY2sKICBjb2xsYXBzZSA9IEZBTFNFCikKIyBTdGFydCB0aGUgZmlndXJlIGNvdW50ZXIKZmlnX2NvdW50IDwtIDAKIyBEZWZpbmUgdGhlIGNhcHRpb25lciBmdW5jdGlvbgpjYXB0aW9uZXIgPC0gZnVuY3Rpb24oY2FwdGlvbikgewogIGZpZ19jb3VudCA8PC0gZmlnX2NvdW50ICsgMQogIHBhc3RlMCgiRmlndXJlICIsIGZpZ19jb3VudCwgIjogIiwgY2FwdGlvbikKfQojIERlZmluZSB0aGUgZnVuY3Rpb24gdG8gdHJ1bmNhdGUgYSBudW1iZXIgdG8gdHdvIGRlY2ltYWwgcGxhY2VzCnRydW5jYXRlX3RvX3R3byA8LSBmdW5jdGlvbih4KSB7CiAgZmxvb3IoeCAqIDEwMCkgLyAxMDAKfQoKbTF0YWJsZSA8LSByU1BERTo6Om0xdGFibGUKbTJ0YWJsZSA8LSByU1BERTo6Om0ydGFibGUKbTN0YWJsZSA8LSByU1BERTo6Om0zdGFibGUKbTR0YWJsZSA8LSByU1BERTo6Om00dGFibGUKYGBgCgoKYGBge3J9CiMgaW5zdGFsbC5wYWNrYWdlcygiSU5MQSIscmVwb3M9YyhnZXRPcHRpb24oInJlcG9zIiksSU5MQT0iaHR0cHM6Ly9pbmxhLnItaW5sYS1kb3dubG9hZC5vcmcvUi90ZXN0aW5nIiksIGRlcD1UUlVFKQojIGlubGEudXBncmFkZSh0ZXN0aW5nID0gVFJVRSkKIyByZW1vdGVzOjppbnN0YWxsX2dpdGh1YigiaW5sYWJydS1vcmcvaW5sYWJydSIsIHJlZiA9ICJkZXZlbCIpCiMgcmVtb3Rlczo6aW5zdGFsbF9naXRodWIoImRhdmlkYm9saW4vcnNwZGUiLCByZWYgPSAiZGV2ZWwiKQojIHJlbW90ZXM6Omluc3RhbGxfZ2l0aHViKCJkYXZpZGJvbGluL21ldHJpY2dyYXBoIiwgcmVmID0gImRldmVsIikKbGlicmFyeShJTkxBKQpsaWJyYXJ5KGlubGFicnUpCmxpYnJhcnkoclNQREUpCmxpYnJhcnkoTWV0cmljR3JhcGgpCmxpYnJhcnkoZ3JhdGVmdWwpCgpsaWJyYXJ5KHBsb3RseSkKYGBgCgoKIyBVdGlsaXRhcnkgZnVuY3Rpb25zCgpgYGB7cn0KIyBGb3IgZWFjaCBtIGFuZCBiZXRhLCB0aGlzIGZ1bmN0aW9uIHJldHVybnMgY19tL2Jfe20rMX0gYW5kIHRoZSByb290cyBvZiByYiBhbmQgcmMKbXkuZ2V0LnJvb3RzIDwtIGZ1bmN0aW9uKG9yZGVyLCBiZXRhKSB7CiAgbXQgPC0gZ2V0KHBhc3RlMCgibSIsIG9yZGVyLCAidGFibGUiKSkKICByYiA8LSByZXAoMCwgb3JkZXIgKyAxKQogIHJjIDwtIHJlcCgwLCBvcmRlcikKICBpZihvcmRlciA9PSAxKSB7CiAgICAgIHJjID0gYXBwcm94KG10JGJldGEsIG10W1twYXN0ZTAoInJjIildXSwgYmV0YSkkeQogICAgfSBlbHNlIHsKICAgICAgcmMgPSBzYXBwbHkoMTpvcmRlciwgZnVuY3Rpb24oaSkgewogICAgICAgIGFwcHJveChtdCRiZXRhLCBtdFtbcGFzdGUwKCJyYy4iLCBpKV1dLCBiZXRhKSR5CiAgICAgIH0pCiAgICB9CiAgICByYiA9IHNhcHBseSgxOihvcmRlcisxKSwgZnVuY3Rpb24oaSkgewogICAgICBhcHByb3gobXQkYmV0YSwgbXRbW3Bhc3RlMCgicmIuIiwgaSldXSwgeG91dCA9IGJldGEpJHkKICAgIH0pCiAgICBmYWN0b3IgPSBhcHByb3gobXQkYmV0YSwgbXQkZmFjdG9yLCB4b3V0ID0gYmV0YSkkeQogIHJldHVybihsaXN0KHJiID0gcmIsIHJjID0gcmMsIGZhY3RvciA9IGZhY3RvcikpCn0KCiMgRnVuY3Rpb24gdGhlIHBvbHlub21pYWwgY29lZmZpY2llbnRzIGluIGluY3JlYXNpbmcgb3JkZXIgbGlrZSBhK2J4K2N4XjIrLi4uCnBvbHkuZnJvbS5yb290cyA8LSBmdW5jdGlvbihyb290cykgewogIAogIGNvZWYgPC0gMQogIGZvciAociBpbiByb290cykge2NvZWYgPC0gY29udm9sdmUoY29lZiwgYygxLCAtciksIHR5cGUgPSAib3BlbiIpfQogIHJldHVybihjb2VmKQp9CiMgRnVuY3Rpb24gdG8gY29tcHV0ZSB0aGUgcGFydGlhbCBmcmFjdGlvbiBwYXJhbWV0ZXJzCmNvbXB1dGUucGFydGlhbC5mcmFjdGlvbi5wYXJhbSA8LSBmdW5jdGlvbihmYWN0b3IsIHByX3Jvb3RzLCBwbF9yb290cywgY3RlKSB7CiAgCiAgcHJfY29lZiA8LSBjKDAsIHBvbHkuZnJvbS5yb290cyhwcl9yb290cykpIAogIHBsX2NvZWYgPC0gcG9seS5mcm9tLnJvb3RzKHBsX3Jvb3RzKSAKICBmYWN0b3JfcHJfY29lZiA8LSBwcl9jb2VmCiAgcHJfcGx1c19wbF9jb2VmIDwtIGZhY3Rvcl9wcl9jb2VmICsgY3RlL2ZhY3RvciAqIHBsX2NvZWYKICByZXMgPC0gZ3NpZ25hbDo6cmVzaWR1ZShmYWN0b3JfcHJfY29lZiwgcHJfcGx1c19wbF9jb2VmKQogIHJldHVybihsaXN0KHIgPSByZXMkciwgcCA9IHJlcyRwLCBrID0gcmVzJGspKSAKfQojIEZ1bmN0aW9uIHRvIGNvbXB1dGUgdGhlIGZyYWN0aW9uYWwgb3BlcmF0b3IKbXkuZnJhY3Rpb25hbC5vcGVyYXRvcnMuZnJhYyA8LSBmdW5jdGlvbihMLCBiZXRhLCBDLCBzY2FsZS5mYWN0b3IsIG0gPSAxLCB0aW1lX3N0ZXApIHsKICAKICBDIDwtIE1hdHJpeDo6RGlhZ29uYWwoZGltKEMpWzFdLCByb3dTdW1zKEMpKSAKICBDaSA8LSBNYXRyaXg6OkRpYWdvbmFsKGRpbShDKVsxXSwgMSAvIHJvd1N1bXMoQykpIAogIEkgPC0gTWF0cml4OjpEaWFnb25hbChkaW0oQylbMV0pCiAgTCA8LSBMIC8gc2NhbGUuZmFjdG9yIAogIExDaSA8LSBMICUqJSBDaQogIGlmKGJldGEgPT0gMSl7CiAgICByZXR1cm4obGlzdChDaSA9IENpLCBDID0gQywgTENpID0gTENpLCBMID0gTCwgbSA9IG0sIGJldGEgPSBiZXRhLCBMSFMgPSBDICsgdGltZV9zdGVwICogTCkpCiAgfSBlbHNlIHsKICAgIHJvb3RzIDwtIG15LmdldC5yb290cyhtLCBiZXRhKQogICAgcG9sZXNfcnNfayA8LSBjb21wdXRlLnBhcnRpYWwuZnJhY3Rpb24ucGFyYW0ocm9vdHMkZmFjdG9yLCByb290cyRyYywgcm9vdHMkcmIsIHRpbWVfc3RlcCkKCiAgICBwYXJ0aWFsX2ZyYWN0aW9uX3Rlcm1zIDwtIGxpc3QoKQogICAgZm9yIChpIGluIDE6KG0rMSkpIHtwYXJ0aWFsX2ZyYWN0aW9uX3Rlcm1zW1tpXV0gPC0gKExDaSAtIHBvbGVzX3JzX2skcFtpXSAqIEkpL3BvbGVzX3JzX2skcltpXX0KICAgIHBhcnRpYWxfZnJhY3Rpb25fdGVybXNbW20rMl1dIDwtIGlmZWxzZShpcy5udWxsKHBvbGVzX3JzX2skayksIDAsIHBvbGVzX3JzX2skaykgKiBJCiAgICByZXR1cm4obGlzdChDaSA9IENpLCBDID0gQywgTENpID0gTENpLCBMID0gTCwgbSA9IG0sIGJldGEgPSBiZXRhLCBwYXJ0aWFsX2ZyYWN0aW9uX3Rlcm1zID0gcGFydGlhbF9mcmFjdGlvbl90ZXJtcykpCiAgfQp9CiMgRnVuY3Rpb24gdG8gc29sdmUgdGhlIGl0ZXJhdGlvbgpteS5zb2x2ZXIuZnJhYyA8LSBmdW5jdGlvbihvYmosIHYpewogIGJldGEgPC0gb2JqJGJldGEKICBtIDwtIG9iaiRtCiAgQyA8LSBvYmokQwogIENpIDwtIG9iaiRDaQogIGlmIChiZXRhID09IDEpewogICAgcmV0dXJuKHNvbHZlKG9iaiRMSFMsIHYpKQogIH0gZWxzZSB7CiAgICBwYXJ0aWFsX2ZyYWN0aW9uX3Rlcm1zIDwtIG9iaiRwYXJ0aWFsX2ZyYWN0aW9uX3Rlcm1zCiAgICBvdXRwdXQgPC0gcGFydGlhbF9mcmFjdGlvbl90ZXJtc1tbbSsyXV0gJSolIHYKICAgIGZvciAoaSBpbiAxOihtKzEpKSB7b3V0cHV0IDwtIG91dHB1dCArIHNvbHZlKHBhcnRpYWxfZnJhY3Rpb25fdGVybXNbW2ldXSwgdil9CiAgICByZXR1cm4oQ2kgJSolIG91dHB1dCkKICB9Cn0KIyBGdW5jdGlvbiB0byBidWlsZCBhIHRhZHBvbGUgZ3JhcGggYW5kIGNyZWF0ZSBhIG1lc2gKZ2V0cy5ncmFwaC50YWRwb2xlIDwtIGZ1bmN0aW9uKGgpewogIGVkZ2UxIDwtIHJiaW5kKGMoMCwwKSxjKDEsMCkpCiAgdGhldGEgPC0gc2VxKGZyb209LXBpLHRvPXBpLGxlbmd0aC5vdXQgPSAxMDAwMCkKICBlZGdlMiA8LSBjYmluZCgxKzEvcGkrY29zKHRoZXRhKS9waSxzaW4odGhldGEpL3BpKQogIGVkZ2VzIDwtIGxpc3QoZWRnZTEsIGVkZ2UyKQogIGdyYXBoIDwtIG1ldHJpY19ncmFwaCRuZXcoZWRnZXMgPSBlZGdlcykKICBncmFwaCRzZXRfbWFudWFsX2VkZ2VfbGVuZ3RocyhlZGdlX2xlbmd0aHMgPSBjKDEsMikpCiAgZ3JhcGgkYnVpbGRfbWVzaChoID0gaCkKICByZXR1cm4oZ3JhcGgpCn0KIyBGdW5jdGlvbiB0byBjb21wdXRlIHRoZSBlaWdlbmZ1bmN0aW9ucyAKdGFkcG9sZS5laWcgPC0gZnVuY3Rpb24oayxncmFwaCl7CngxIDwtIGMoMCxncmFwaCRnZXRfZWRnZV9sZW5ndGhzKClbMV0qZ3JhcGgkbWVzaCRQdEVbZ3JhcGgkbWVzaCRQdEVbLDFdPT0xLDJdKSAKeDIgPC0gYygwLGdyYXBoJGdldF9lZGdlX2xlbmd0aHMoKVsyXSpncmFwaCRtZXNoJFB0RVtncmFwaCRtZXNoJFB0RVssMV09PTIsMl0pIAoKaWYoaz09MCl7IAogIGYuZTEgPC0gcmVwKDEsbGVuZ3RoKHgxKSkgCiAgZi5lMiA8LSByZXAoMSxsZW5ndGgoeDIpKSAKICBmMSA9IGMoZi5lMVsxXSxmLmUyWzFdLGYuZTFbLTFdLCBmLmUyWy0xXSkgCiAgZiA9IGxpc3QocGhpPWYxL3NxcnQoMykpIAogIAp9IGVsc2UgewogIGYuZTEgPC0gLTIqc2luKHBpKmsqMS8yKSpjb3MocGkqayp4MS8yKSAKICBmLmUyIDwtIHNpbihwaSprKngyLzIpICAgICAgICAgICAgICAgICAgCiAgCiAgZjEgPSBjKGYuZTFbMV0sZi5lMlsxXSxmLmUxWy0xXSwgZi5lMlstMV0pIAogIAogIGlmKChrICUlIDIpPT0xKXsgCiAgICBmID0gbGlzdChwaGk9ZjEvc3FydCgzKSkgCiAgfSBlbHNlIHsgCiAgICBmLmUxIDwtICgtMSlee2svMn0qY29zKHBpKmsqeDEvMikKICAgIGYuZTIgPC0gY29zKHBpKmsqeDIvMikKICAgIGYyID0gYyhmLmUxWzFdLGYuZTJbMV0sZi5lMVstMV0sZi5lMlstMV0pIAogICAgZiA8LSBsaXN0KHBoaT1mMSxwc2k9ZjIvc3FydCgzLzIpKQogIH0KfQpyZXR1cm4oZikKfQojIEZ1bmN0aW9uIHRvIG9yZGVyIHRoZSB2ZXJ0aWNlcyBmb3IgcGxvdHRpbmcKcGxvdHRpbmcub3JkZXIgPC0gZnVuY3Rpb24odiwgZ3JhcGgpewogIGVkZ2VfbnVtYmVyIDwtIGdyYXBoJG1lc2gkVnRFWywgMV0KICBwb3MgPC0gc3VtKGVkZ2VfbnVtYmVyID09IDEpKzEKICByZXR1cm4oYyh2WzFdLCB2WzM6cG9zXSwgdlsyXSwgdlsocG9zKzEpOmxlbmd0aCh2KV0sIHZbMl0pKQp9CiMgRnVuY3Rpb24gdG8gcGxvdCBpbiAzRApncmFwaC5wbG90dGVyLjNkIDwtIGZ1bmN0aW9uKGdyYXBoLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICBVX3RydWUsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIFVfYXBwcm94MSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgI1VfYXBwcm94MiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdGltZV9zZXEsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIHRpbWVfc3RlcCl7CiAgeCA8LSBncmFwaCRtZXNoJFZbLCAxXQogIHkgPC0gZ3JhcGgkbWVzaCRWWywgMl0KICB4IDwtIHBsb3R0aW5nLm9yZGVyKHgsIGdyYXBoKQogIHkgPC0gcGxvdHRpbmcub3JkZXIoeSwgZ3JhcGgpCiAgd2VpZ2h0cyA8LSBncmFwaCRtZXNoJHdlaWdodHMKICAKICBjdW1zdW0xIDwtIHNxcnQodGltZV9zdGVwICogY3Vtc3VtKHQod2VpZ2h0cykgJSolIChVX3RydWUgLSBVX2FwcHJveDEpXjIpKQogICNjdW1zdW0yIDwtIHNxcnQodGltZV9zdGVwICogY3Vtc3VtKHQod2VpZ2h0cykgJSolKFVfdHJ1ZSAtIFVfYXBwcm94MileMikpCiAgI2N1bXN1bTMgPC0gc3FydCh0aW1lX3N0ZXAgKiBjdW1zdW0odCh3ZWlnaHRzKSAlKiUgKFVfYXBwcm94MSAtIFVfYXBwcm94MileMikpCiAgCiAgVV90cnVlIDwtIGFwcGx5KFVfdHJ1ZSwgMiwgcGxvdHRpbmcub3JkZXIsIGdyYXBoID0gZ3JhcGgpCiAgVV9hcHByb3gxIDwtIGFwcGx5KFVfYXBwcm94MSwgMiwgcGxvdHRpbmcub3JkZXIsIGdyYXBoID0gZ3JhcGgpCiAgI1VfYXBwcm94MiA8LSBhcHBseShVX2FwcHJveDIsIDIsIHBsb3R0aW5nLm9yZGVyLCBncmFwaCA9IGdyYXBoKQogIAogIHBsb3RfZGF0YSA8LSBkYXRhLmZyYW1lKAogIHggPSByZXAoeCwgdGltZXMgPSBuY29sKFVfdHJ1ZSkpLAogIHkgPSByZXAoeSwgdGltZXMgPSBuY29sKFVfdHJ1ZSkpLAogIHpfdHJ1ZSA9IGFzLnZlY3RvcihVX3RydWUpLAogIHpfYXBwcm94MSA9IGFzLnZlY3RvcihVX2FwcHJveDEpLAogICN6X2FwcHJveDIgPSBhcy52ZWN0b3IoVV9hcHByb3gyKSwKICBmcmFtZSA9IHJlcCh0aW1lX3NlcSwgZWFjaCA9IGxlbmd0aCh4KSkpCiAgCiAgeF9yYW5nZSA8LSByYW5nZSh4KQogIHlfcmFuZ2UgPC0gcmFuZ2UoeSkKICB6X3JhbmdlIDwtIHJhbmdlKGMoVV90cnVlLCAKICAgICAgICAgICAgICAgICAgICAgI1VfYXBwcm94MiwgCiAgICAgICAgICAgICAgICAgICAgIFVfYXBwcm94MSkpCiAgCiAgcCA8LSBwbG90X2x5KHBsb3RfZGF0YSwgZnJhbWUgPSB+ZnJhbWUpICU+JQogIGFkZF90cmFjZSh4ID0gfngsIHkgPSB+eSwgeiA9IH56X3RydWUsIHR5cGUgPSAic2NhdHRlcjNkIiwgbW9kZSA9ICJsaW5lcyIsIG5hbWUgPSAiVCIsCiAgICAgICAgICAgIGxpbmUgPSBsaXN0KGNvbG9yID0gImdyZWVuIiwgd2lkdGggPSAyKSkgJT4lCiAgYWRkX3RyYWNlKHggPSB+eCwgeSA9IH55LCB6ID0gfnpfYXBwcm94MSwgdHlwZSA9ICJzY2F0dGVyM2QiLCBtb2RlID0gImxpbmVzIiwgbmFtZSA9ICJBMSIsCiAgICAgICAgICAgIGxpbmUgPSBsaXN0KGNvbG9yID0gInJlZCIsIHdpZHRoID0gMikpICU+JQogICMgYWRkX3RyYWNlKHggPSB+eCwgeSA9IH55LCB6ID0gfnpfYXBwcm94MiwgdHlwZSA9ICJzY2F0dGVyM2QiLCBtb2RlID0gImxpbmVzIiwgbmFtZSA9ICJBMiIsIAogICMgICAgICAgICAgIGxpbmUgPSBsaXN0KGNvbG9yID0gImJsdWUiLCB3aWR0aCA9IDIpKSAlPiUKICBsYXlvdXQoCiAgICBzY2VuZSA9IGxpc3QoCiAgICAgIHhheGlzID0gbGlzdCh0aXRsZSA9ICJ4IiwgcmFuZ2UgPSB4X3JhbmdlKSwKICAgICAgeWF4aXMgPSBsaXN0KHRpdGxlID0gInkiLCByYW5nZSA9IHlfcmFuZ2UpLAogICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAieiIsIHJhbmdlID0gel9yYW5nZSksCiAgICAgIGFzcGVjdHJhdGlvID0gbGlzdCh4ID0gMi40LCB5ID0gMS4yLCB6ID0gMS4yKSwKICAgICAgY2FtZXJhID0gbGlzdChleWUgPSBsaXN0KHggPSAxLjUsIHkgPSAxLjUsIHogPSAxKSwgY2VudGVyID0gbGlzdCh4ID0gMCwgeSA9IDAsIHogPSAwKSkpLAogICAgdXBkYXRlbWVudXMgPSBsaXN0KGxpc3QodHlwZSA9ICJidXR0b25zIiwgc2hvd2FjdGl2ZSA9IEZBTFNFLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgYnV0dG9ucyA9IGxpc3QoCiAgICAgICAgICBsaXN0KGxhYmVsID0gIlBsYXkiLCBtZXRob2QgPSAiYW5pbWF0ZSIsCiAgICAgICAgICAgICAgIGFyZ3MgPSBsaXN0KE5VTEwsIGxpc3QoZnJhbWUgPSBsaXN0KGR1cmF0aW9uID0gMTAwLCByZWRyYXcgPSBUUlVFKSwgZnJvbWN1cnJlbnQgPSBUUlVFKSkpLAogICAgICAgICAgbGlzdChsYWJlbCA9ICJQYXVzZSIsIG1ldGhvZCA9ICJhbmltYXRlIiwKICAgICAgICAgICAgICAgYXJncyA9IGxpc3QoTlVMTCwgbGlzdChtb2RlID0gImltbWVkaWF0ZSIsIGZyYW1lID0gbGlzdChkdXJhdGlvbiA9IDApLCByZWRyYXcgPSBGQUxTRSkpKQogICAgICAgICkKICAgICAgKQogICAgKSwKICAgIHRpdGxlID0gIlRpbWU6IDAiCiAgKSAlPiUgcGxvdGx5X2J1aWxkKCkKICBmb3IgKGkgaW4gc2VxX2Fsb25nKHAkeCRmcmFtZXMpKSB7CiAgICB0IDwtIHRpbWVfc2VxW2ldCiAgICBlcnIxIDwtIHNpZ25pZihjdW1zdW0xW2ldLCA2KQogICAgI2VycjIgPC0gc2lnbmlmKGN1bXN1bTJbaV0sIDQpCiAgICAjZXJyMyA8LSBzaWduaWYoY3Vtc3VtM1tpXSwgNCkKICAgIHAkeCRmcmFtZXNbW2ldXSRsYXlvdXQgPC0gbGlzdCh0aXRsZSA9IHBhc3RlMCgiVGltZTogIiwgdCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyAiIHwgY3VtIEU9VC1BMjogIiwgZXJyMiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyAiIHwgY3VtIEU9QTEtQTI6ICIsIGVycjMsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIiB8IGN1bSBFPVQtQTE6ICIsIGVycjEpKQogIH0KICByZXR1cm4ocCkKfQojIEZ1bmN0aW9uIHRvIHBsb3QgdGhlIGVycm9yIGF0IGVhY2ggdGltZSBzdGVwCmVycm9yLmF0LmVhY2gudGltZS5wbG90dGVyIDwtIGZ1bmN0aW9uKGdyYXBoLCBVX3RydWUsIFVfYXBwcm94MSwgVV9hcHByb3gyLCB0aW1lX3NlcSwgdGltZV9zdGVwKSB7CiAgd2VpZ2h0cyA8LSBncmFwaCRtZXNoJHdlaWdodHMKICBlcnJvcl9hdF9lYWNoX3RpbWUxIDwtIHQod2VpZ2h0cykgJSolIChVX3RydWUgLSBVX2FwcHJveDEpXjIKICBlcnJvcl9hdF9lYWNoX3RpbWUyIDwtIHQod2VpZ2h0cykgJSolIChVX3RydWUgLSBVX2FwcHJveDIpXjIKICBlcnJvcl9iZXR3ZWVuX2JvdGhfYXBwcm94IDwtIHQod2VpZ2h0cykgJSolIChVX2FwcHJveDEgLSBVX2FwcHJveDIpXjIKICAKICBlcnJvcjEgPC0gc3FydChhcy5kb3VibGUodCh3ZWlnaHRzKSAlKiUgKFVfdHJ1ZSAtIFVfYXBwcm94MSleMiAlKiUgcmVwKHRpbWVfc3RlcCwgbmNvbChVX3RydWUpKSkpCiAgZXJyb3IyIDwtIHNxcnQoYXMuZG91YmxlKHQod2VpZ2h0cykgJSolIChVX3RydWUgLSBVX2FwcHJveDIpXjIgJSolIHJlcCh0aW1lX3N0ZXAsIG5jb2woVV90cnVlKSkpKQogIGVycm9yYiA8LSBzcXJ0KGFzLmRvdWJsZSh0KHdlaWdodHMpICUqJSAoVV9hcHByb3gxIC0gVV9hcHByb3gyKV4yICUqJSByZXAodGltZV9zdGVwLCBuY29sKFVfdHJ1ZSkpKSkKICAKICBmaWcgPC0gcGxvdF9seSgpICU+JSAKICBhZGRfdHJhY2UoCiAgeCA9IH50aW1lX3NlcSwgeSA9IH5lcnJvcl9hdF9lYWNoX3RpbWUxLCB0eXBlID0gJ3NjYXR0ZXInLCBtb2RlID0gJ2xpbmVzK21hcmtlcnMnLAogIGxpbmUgPSBsaXN0KGNvbG9yID0gJ3JlZCcsIHdpZHRoID0gMiksCiAgbWFya2VyID0gbGlzdChjb2xvciA9ICdyZWQnLCBzaXplID0gNCksCiAgbmFtZSA9IHBhc3RlMCgiRT1ULUExOiAiLCBzcHJpbnRmKCIlLjNlIiwgZXJyb3IxKSkKKSAlPiUgCiAgYWRkX3RyYWNlKAogIHggPSB+dGltZV9zZXEsIHkgPSB+ZXJyb3JfYXRfZWFjaF90aW1lMiwgdHlwZSA9ICdzY2F0dGVyJywgbW9kZSA9ICdsaW5lcyttYXJrZXJzJywKICBsaW5lID0gbGlzdChjb2xvciA9ICdibHVlJywgd2lkdGggPSAyLCBkYXNoID0gImRvdCIpLAogIG1hcmtlciA9IGxpc3QoY29sb3IgPSAnYmx1ZScsIHNpemUgPSA0KSwKICBuYW1lID0gcGFzdGUwKCJFPVQtQTI6ICIsIHNwcmludGYoIiUuM2UiLCBlcnJvcjIpKQopICU+JSAKICBhZGRfdHJhY2UoCiAgeCA9IH50aW1lX3NlcSwgeSA9IH5lcnJvcl9iZXR3ZWVuX2JvdGhfYXBwcm94LCB0eXBlID0gJ3NjYXR0ZXInLCBtb2RlID0gJ2xpbmVzK21hcmtlcnMnLAogIGxpbmUgPSBsaXN0KGNvbG9yID0gJ29yYW5nZScsIHdpZHRoID0gMiksCiAgbWFya2VyID0gbGlzdChjb2xvciA9ICdvcmFuZ2UnLCBzaXplID0gNCksCiAgbmFtZSA9IHBhc3RlMCgiRT1BMS1BMjogIiwgc3ByaW50ZigiJS4zZSIsIGVycm9yYikpCikgJT4lIAogIGxheW91dCgKICB0aXRsZSA9ICJFcnJvciBhdCBFYWNoIFRpbWUgU3RlcCIsCiAgeGF4aXMgPSBsaXN0KHRpdGxlID0gIkVycm9yIGF0IEVhY2ggVGltZSBTdGVwIiksCiAgeWF4aXMgPSBsaXN0KHRpdGxlID0gIkVycm9yIiksCiAgbGVnZW5kID0gbGlzdCh4ID0gMC4xLCB5ID0gMC45KQopCiAgcmV0dXJuKGZpZykKfQojIEZ1bmN0aW9uIHRvIGNvbXB1dGUgdGhlIGVpZ2VudmFsdWVzIGFuZCBlaWdlbmZ1bmN0aW9ucwpnZXRzLmVpZ2VuLnBhcmFtcyA8LSBmdW5jdGlvbihOX2Zpbml0ZSA9IDQsIGthcHBhID0gMSwgYWxwaGEgPSAwLjUsIGdyYXBoKXsKICBFSUdFTlZBTF9BTFBIQSA8LSBOVUxMCiAgRUlHRU5GVU4gPC0gTlVMTAogIGZvciAoaiBpbiAwOk5fZmluaXRlKSB7CiAgICAgIGxhbWJkYV9qX2FscGhhIDwtIChrYXBwYV4yICsgKGoqcGkvMileMileKGFscGhhLzIpCiAgICAgIGVfaiA8LSB0YWRwb2xlLmVpZyhqLGdyYXBoKSRwaGkKICAgICAgRUlHRU5WQUxfQUxQSEEgPC0gYyhFSUdFTlZBTF9BTFBIQSwgbGFtYmRhX2pfYWxwaGEpICAgICAgICAgCiAgICAgIEVJR0VORlVOIDwtIGNiaW5kKEVJR0VORlVOLCBlX2opCiAgICAgIGlmIChqPjAgJiYgKGogJSUgMiA9PSAwKSkgewogICAgICAgIGxhbWJkYV9qX2FscGhhIDwtIChrYXBwYV4yICsgKGoqcGkvMileMileKGFscGhhLzIpCiAgICAgICAgZV9qIDwtIHRhZHBvbGUuZWlnKGosZ3JhcGgpJHBzaQogICAgICAgIEVJR0VOVkFMX0FMUEhBIDwtIGMoRUlHRU5WQUxfQUxQSEEsIGxhbWJkYV9qX2FscGhhKSAgICAgICAgIAogICAgICAgIEVJR0VORlVOIDwtIGNiaW5kKEVJR0VORlVOLCBlX2opCiAgICAgIH0KICB9CiAgcmV0dXJuKGxpc3QoRUlHRU5WQUxfQUxQSEEgPSBFSUdFTlZBTF9BTFBIQSwKICAgICAgICAgICAgICBFSUdFTkZVTiA9IEVJR0VORlVOKSkKfQpgYGAKCgpXZSB3YW50IHRvIHNvbHZlIHRoZSBmcmFjdGlvbmFsIGRpZmZ1c2lvbiBlcXVhdGlvbgpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTptYWluZXF9CiAgICBccGFydGlhbF90IHUrKFxrYXBwYV4yLVxEZWx0YV9cR2FtbWEpXntcYWxwaGEvMn0gdT1mIFx0ZXh0IHsgb24gfSBcR2FtbWEgXHRpbWVzKDAsIFQpLCBccXVhZCB1KDApPXVfMCBcdGV4dCB7IG9uIH0gXEdhbW1hLApcZW5ke2VxdWF0aW9ufQp3aGVyZSAkdSQgc2F0aXNmaWVzIHRoZSBLaXJjaGhvZmYgdmVydGV4IGNvbmRpdGlvbnMKXGJlZ2lue2VxdWF0aW9ufQpcbGFiZWx7ZXE6S2NvbmR9CiAgICBcbGVmdFx7XHBoaVxpbiBDKFxHYW1tYSlcO1xCaWd8XDsgXGZvcmFsbCB2XGluIFY6IFxzdW1fe2VcaW5cbWF0aGNhbHtFfV92fVxwYXJ0aWFsX2UgXHBoaSh2KT0wIFxyaWdodFx9ClxlbmR7ZXF1YXRpb259ClRoZSBzb2x1dGlvbiBpcyBnaXZlbiBieQpcYmVnaW57ZXF1YXRpb259ClxsYWJlbHtlcTpzb2xfcmVwcmVudGF0aW9ufQogICAgICAgIHUocyx0KSA9IFxkaXNwbGF5c3R5bGVcc3VtX3tqXGluXG1hdGhiYntOfX1lXnstXGxhbWJkYV57XGFscGhhLzJ9X2p0fVxsZWZ0KHVfMCwgZV9qXHJpZ2h0KV97TF8yKFxHYW1tYSl9ZV9qKHMpICsgXGludF8wXnQgXGRpc3BsYXlzdHlsZVxzdW1fe2pcaW5cbWF0aGJie059fWVeey1cbGFtYmRhXntcYWxwaGEvMn1faih0LXIpfVxsZWZ0KGYoXGNkb3QsIHIpLCBlX2pccmlnaHQpX3tMXzIoXEdhbW1hKX1lX2oocylkci4KXGVuZHtlcXVhdGlvbn0KCklmIHdlIGNob29zZSAkd19qJCBhbmQgJHZfaiQgYW5kIHRha2UgdGhlIGluaXRpYWwgY29uZGl0aW9uIGFuZCB0aGUgcmlnaHQgaGFuZCBzaWRlIGZ1bmNpdG9uIGFzIAoKXGJlZ2lue2VxdWF0aW9ufQogICAgdV8wKHMpID0gXHN1bV97aj0wfV57Tn0gd19qIGVfaihzKSBcdGV4dHsgYW5kIHNvIH0gXGxlZnQodV8wLCBlX2pccmlnaHQpX3tMXzIoXEdhbW1hKX0gPSB3X2osClxlbmR7ZXF1YXRpb259ClxiZWdpbntlcXVhdGlvbn0KICAgXHRleHR7SW4gbWF0cml4IG5vdGF0aW9uOiB9IFxxdWFkXGJvbGRzeW1ib2x7VX1fMCA9IFxib2xkc3ltYm9se0V9X2hcYm9sZHN5bWJvbHtjfSwgXHF1YWQgIFxib2xkc3ltYm9se0V9X2ggPSBcbGVmdFtlXzAsIGVfMSwgXGxkb3RzLCBlX3tOfVxyaWdodF0sIFxxdWFkIFxib2xkc3ltYm9se3d9ID0gXGxlZnRbd18wLCB3XzEsIFxsZG90cywgd197Tn1ccmlnaHRdXlx0b3AsClxlbmR7ZXF1YXRpb259ClxiZWdpbntlcXVhdGlvbn0KICAgXHRleHR7SW4gfSBcdGV4dHR0e1J9OiBcdGV4dHR0e1VfMCA8LSBFSUdFTkZVTiAlKiUgY29lZmZfVV8wfQpcZW5ke2VxdWF0aW9ufQphbmQKXGJlZ2lue2VxdWF0aW9ufQogICAgZihzLHQpID0gXHN1bV97aj0wfV57Tn0gdl9qIGVeey1cbGFtYmRhXntcYWxwaGEvMn1fanR9IGVfaihzKSBcdGV4dHsgYW5kIHNvIH0gXGxlZnQoZihcY2RvdCxyKSwgZV9qXHJpZ2h0KV97TF8yKFxHYW1tYSl9ID0gdl9qIGVeey1cbGFtYmRhXntcYWxwaGEvMn1fanJ9LApcZW5ke2VxdWF0aW9ufQpcYmVnaW57ZXF1YXRpb259CiAgIFx0ZXh0e0luIG1hdHJpeCBub3RhdGlvbjogfSBccXVhZFxib2xkc3ltYm9se2Z9ID0gXGJvbGRzeW1ib2x7RX1faCBcYm9sZHN5bWJvbHtWfSwgXHF1YWQgXGJvbGRzeW1ib2x7Vn1fe2ppfSA9IHZfamVeey1cbGFtYmRhXntcYWxwaGEvMn1fanRfaX0KXGVuZHtlcXVhdGlvbn0KXGJlZ2lue2VxdWF0aW9ufQogICBcdGV4dHtJbiB9IFx0ZXh0dHR7Un06IFx0ZXh0dHR7RkYgPC0gRUlHRU5GVU4gJSolIChjb2VmZl9GRiAqIGV4cCgtb3V0ZXIoRUlHRU5WQUxfQUxQSEEsIHRpbWVfc2VxKSkpfQpcZW5ke2VxdWF0aW9ufQp0aGVuIHRoZSBzb2x1dGlvbiBpcyBnaXZlbiBieQpcYmVnaW57YWxpZ259CiAgICAgICAgdShzLHQpICY9IFxkaXNwbGF5c3R5bGVcc3VtX3tqPTB9XntOfXdfamVeey1cbGFtYmRhXntcYWxwaGEvMn1fanR9ZV9qKHMpICsgXGludF8wXnQgXGRpc3BsYXlzdHlsZVxzdW1fe2o9MH1eTmVeey1cbGFtYmRhXntcYWxwaGEvMn1faih0LXIpfXZfaiBlXnstXGxhbWJkYV57XGFscGhhLzJ9X2pyfWVfaihzKWRyXFwKICAgICAgICAmPSBcZGlzcGxheXN0eWxlXHN1bV97aj0wfV57Tn13X2plXnstXGxhbWJkYV57XGFscGhhLzJ9X2p0fWVfaihzKSArIHQgXGRpc3BsYXlzdHlsZVxzdW1fe2o9MH1eTnZfaiBlXnstXGxhbWJkYV57XGFscGhhLzJ9X2p0fWVfaihzKVxcCiAgICAgICAgJj0gXGRpc3BsYXlzdHlsZVxzdW1fe2o9MH1ee059d19qZV57LVxsYW1iZGFee1xhbHBoYS8yfV9qdH1lX2oocykgKyB0ZihzLHQpXFwKICAgICAgICAmPSBcZGlzcGxheXN0eWxlXHN1bV97aj0wfV57Tn0od19qK3R2X2opZV57LVxsYW1iZGFee1xhbHBoYS8yfV9qdH1lX2oocykgClxlbmR7YWxpZ259CgpcYmVnaW57ZXF1YXRpb259CiAgIFx0ZXh0e0luIG1hdHJpeCBub3RhdGlvbjogfSBccXVhZFxib2xkc3ltYm9se1V9ID1cYm9sZHN5bWJvbHtFfV9oIFxib2xkc3ltYm9se1d9ICsgKFxib2xkc3ltYm9se3R9XGJvbGRzeW1ib2x7Zn1eXHRvcCleXHRvcCwgXHF1YWQgXGJvbGRzeW1ib2x7V31fe2ppfSA9IHdfamVeey1cbGFtYmRhXntcYWxwaGEvMn1fanRfaX0sXHF1YWQgXGJvbGRzeW1ib2x7dH0gPSBcbGVmdFt0XzAsIHRfMSwgXGxkb3RzLCB0X3tLfVxyaWdodF1eXHRvcApcZW5ke2VxdWF0aW9ufQpcYmVnaW57ZXF1YXRpb259CiAgIFx0ZXh0e0luIH0gXHRleHR0dHtSfTogXHRleHR0dHtVX3RydWUgPC0gRUlHRU5GVU4gJSolIChjb2VmZl9VXzAgKiBleHAoLW91dGVyKEVJR0VOVkFMX0FMUEhBLCB0aW1lX3NlcSkpKSsgdCh0aW1lX3NlcSAqIHQoRkYpKX0KXGVuZHtlcXVhdGlvbn0KCmBgYHtyfQprYXBwYSA8LSAxCmFscGhhIDwtIDAuNSAjIGZyb20gMC41IHRvIDIKbSA9IDEKYmV0YSA8LSBhbHBoYS8yCk5fZmluaXRlID0gNCAjIGNob29zZSBldmVuCmFkanVzdGVkX05fZmluaXRlIDwtIE5fZmluaXRlICsgTl9maW5pdGUvMiArIDEKIyBDb2VmZmljaWVudHMgZm9yIHVfMCBhbmQgZgpjb2VmZl9VXzAgPC0gNTAqKDE6YWRqdXN0ZWRfTl9maW5pdGUpXi0xCmNvZWZmX1VfMFstNV0gPC0gMApjb2VmZl9GRiA8LSByZXAoMCwgYWRqdXN0ZWRfTl9maW5pdGUpCmNvZWZmX0ZGWzddIDwtIDEwClRfZmluYWwgPC0gMgoKb3ZlcmtpbGxfdGltZV9wb3dlciA8LSA5Cm92ZXJraWxsX2hfcG93ZXIgPC0gOQoKdGltZV9zdGVwcyA8LSAwLjEgKiAyXi1jKG92ZXJraWxsX3RpbWVfcG93ZXI6MCkKaHMgPC0gMC4xICogMl4tYyhvdmVya2lsbF9oX3Bvd2VyOjApCgpvdmVya2lsbF90aW1lX3N0ZXAgPC0gdGltZV9zdGVwc1sxXQpvdmVya2lsbF9oIDwtIGhzWzFdCgpvdmVya2lsbF90aW1lX3NlcSA8LSBzZXEoMCwgVF9maW5hbCwgYnkgPSBvdmVya2lsbF90aW1lX3N0ZXApCgpvdmVya2lsbF9ncmFwaCA8LSBnZXRzLmdyYXBoLnRhZHBvbGUoaCA9IG92ZXJraWxsX2gpCm92ZXJraWxsX2dyYXBoJGNvbXB1dGVfZmVtKCkKb3ZlcmtpbGxfd2VpZ2h0cyA8LSBvdmVya2lsbF9ncmFwaCRtZXNoJHdlaWdodHMKb3ZlcmtpbGxfZWlnZW5fcGFyYW1zIDwtIGdldHMuZWlnZW4ucGFyYW1zKE5fZmluaXRlID0gTl9maW5pdGUsIGthcHBhID0ga2FwcGEsIGFscGhhID0gYWxwaGEsIGdyYXBoID0gb3ZlcmtpbGxfZ3JhcGgpCkVJR0VOVkFMX0FMUEhBIDwtIG92ZXJraWxsX2VpZ2VuX3BhcmFtcyRFSUdFTlZBTF9BTFBIQQpvdmVya2lsbF9FSUdFTkZVTiA8LSBvdmVya2lsbF9laWdlbl9wYXJhbXMkRUlHRU5GVU4KCgpvdmVya2lsbF9VX3RydWUgPC0gb3ZlcmtpbGxfRUlHRU5GVU4gJSolIAogIG91dGVyKDE6bGVuZ3RoKGNvZWZmX1VfMCksIDE6bGVuZ3RoKG92ZXJraWxsX3RpbWVfc2VxKSwgCiAgICAgICAgZnVuY3Rpb24oaSwgaikgKGNvZWZmX1VfMFtpXSArIGNvZWZmX0ZGW2ldICogb3ZlcmtpbGxfdGltZV9zZXFbal0pICogZXhwKC1FSUdFTlZBTF9BTFBIQVtpXSAqIG92ZXJraWxsX3RpbWVfc2VxW2pdKSkKCmJ5X3ZlY3RvciA8LSAyXmMoMDpvdmVya2lsbF90aW1lX3Bvd2VyKQoKZXJyb3JzIDwtIG1hdHJpeChOQSwgbnJvdyA9IGxlbmd0aChocyksIG5jb2wgPSBsZW5ndGgodGltZV9zdGVwcykpCmZvciAoaSBpbiAxOmxlbmd0aChocykpIHsKICBoIDwtIGhzW2ldCiAgZ3JhcGggPC0gZ2V0cy5ncmFwaC50YWRwb2xlKGggPSBoKQogIGdyYXBoJGNvbXB1dGVfZmVtKCkKICBHIDwtIGdyYXBoJG1lc2gkRwogIEMgPC0gZ3JhcGgkbWVzaCRDCiAgTCA8LSBrYXBwYV4yKkMgKyBHCiAgZWlnZW5fcGFyYW1zIDwtIGdldHMuZWlnZW4ucGFyYW1zKE5fZmluaXRlID0gTl9maW5pdGUsIGthcHBhID0ga2FwcGEsIGFscGhhID0gYWxwaGEsIGdyYXBoID0gZ3JhcGgpCiAgRUlHRU5GVU4gPC0gZWlnZW5fcGFyYW1zJEVJR0VORlVOCiAgVV8wIDwtIEVJR0VORlVOICUqJSBjb2VmZl9VXzAKICBBIDwtIGdyYXBoJGZlbV9iYXNpcyhvdmVya2lsbF9ncmFwaCRnZXRfbWVzaF9sb2NhdGlvbnMoKSkKICBmb3IgKGogaW4gMTpsZW5ndGgodGltZV9zdGVwcykpIHsKICAgIHRpbWVfc3RlcCA8LSB0aW1lX3N0ZXBzW2pdCiAgICBwcmludChieV92ZWN0b3Jbal0pCiAgICBjb2Fyc2VfaW5kaWNlcyA8LSBzZXEoMSwgbGVuZ3RoKG92ZXJraWxsX3RpbWVfc2VxKSwgYnkgPSBieV92ZWN0b3Jbal0pCiAgICB0aW1lX3NlcSA8LSBvdmVya2lsbF90aW1lX3NlcVtjb2Fyc2VfaW5kaWNlc10KICAgIG15X29wX2ZyYWMgPC0gbXkuZnJhY3Rpb25hbC5vcGVyYXRvcnMuZnJhYyhMLCBiZXRhLCBDLCBzY2FsZS5mYWN0b3IgPSBrYXBwYV4yLCBtID0gbSwgdGltZV9zdGVwKQogICAgSU5UX0JBU0lTX0VJR0VOIDwtIHQob3ZlcmtpbGxfRUlHRU5GVU4pICUqJSBvdmVya2lsbF9ncmFwaCRtZXNoJEMgJSolIEEKICAgIEZGX2FwcHJveCA8LSB0KElOVF9CQVNJU19FSUdFTikgJSolIAogICAgICBvdXRlcigxOmxlbmd0aChjb2VmZl9GRiksIDE6bGVuZ3RoKHRpbWVfc2VxKSwgCiAgICAgICAgZnVuY3Rpb24oaSwgaikgY29lZmZfRkZbaV0gKiBleHAoLUVJR0VOVkFMX0FMUEhBW2ldICogdGltZV9zZXFbal0pKQogICAgVV9hcHByb3ggPC0gbWF0cml4KE5BLCBucm93ID0gbnJvdyhDKSwgbmNvbCA9IGxlbmd0aCh0aW1lX3NlcSkpCiAgICBVX2FwcHJveFssIDFdIDwtIFVfMAogICAgZm9yIChrIGluIDE6KGxlbmd0aCh0aW1lX3NlcSkgLSAxKSkgewogICAgICBVX2FwcHJveFssIGsgKyAxXSA8LSBhcy5tYXRyaXgobXkuc29sdmVyLmZyYWMobXlfb3BfZnJhYywgbXlfb3BfZnJhYyRDICUqJSBVX2FwcHJveFssIGtdICsgdGltZV9zdGVwICogRkZfYXBwcm94WywgayArIDFdKSkKICAgICAgfQogICAgc2xpY2VkX292ZXJraWxsX1VfdHJ1ZSA8LSBvdmVya2lsbF9VX3RydWVbLCBjb2Fyc2VfaW5kaWNlc10KICAgIHByb2plY3RlZF9VX2FwcHJveCA8LSBBICUqJSBVX2FwcHJveAogICAgZXJyb3JzW2ksal0gPC0gc3FydChhcy5kb3VibGUodChvdmVya2lsbF93ZWlnaHRzKSAlKiUgKHNsaWNlZF9vdmVya2lsbF9VX3RydWUgLSBwcm9qZWN0ZWRfVV9hcHByb3gpXjIgJSolIHJlcCh0aW1lX3N0ZXAsIGxlbmd0aCh0aW1lX3NlcSkpKSkKICB9Cn0KYGBgCgoKYGBge3J9CnByaW50KGVycm9ycykKZ3JhcGgucGxvdHRlci4zZChvdmVya2lsbF9ncmFwaCwgc2xpY2VkX292ZXJraWxsX1VfdHJ1ZSwgcHJvamVjdGVkX1VfYXBwcm94LCB0aW1lX3NlcSwgdGltZV9zdGVwKQpgYGAKCmBgYHtyfQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkocmVzaGFwZTIpCmxpYnJhcnkocGxvdGx5KQoKCiMgQ3JlYXRlIG1hdHJpeAptYXQgPC0gZXJyb3JzCgojIENvbnZlcnQgdG8gZGF0YSBmcmFtZQpkZiA8LSBtZWx0KG1hdCkKY29sbmFtZXMoZGYpIDwtIGMoIlkiLCAiWCIsICJWYWx1ZSIpCgojIGdncGxvdApwIDwtIGdncGxvdChkZiwgYWVzKHggPSBYLCB5ID0gWSwgZmlsbCA9IFZhbHVlLCB0ZXh0ID0gcGFzdGUoIlZhbHVlOiIsIHJvdW5kKFZhbHVlLCAyKSkpKSArCiAgZ2VvbV90aWxlKCkgKwogIHNjYWxlX2ZpbGxfZ3JhZGllbnQobG93ID0gIndoaXRlIiwgaGlnaCA9ICJyZWQiKSArCiAgY29vcmRfZml4ZWQoKSArCiAgdGhlbWVfbWluaW1hbCgpCgpwCgojIENvbnZlcnQgdG8gaW50ZXJhY3RpdmUgcGxvdApnZ3Bsb3RseShwLCB0b29sdGlwID0gInRleHQiKQoKaW1hZ2UoMToxMCwgMToxMCwgdChtYXQpWywgbnJvdyhtYXQpOjFdLCBjb2wgPSBoZWF0LmNvbG9ycygxMDApLCB4bGFiID0gIiIsIHlsYWIgPSAiIikKCmBgYAoK